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Abstract 

The respective roles of local and nonlocal interactions in the thermody- 
namic cooperativity of proteins are investigated using continuum (off-lattice) 
native-centric Go-like models with a coarse-grained C Q chain representation. 
We study a series of models in which the (local) bond- and torsion-angle 
terms have different strengths relative to the (nonlocal) pairwise contact en- 
ergy terms. Conformational distributions in these models are sampled by 
Langevin dynamics. Thermodynamic cooperativity is characterized by the 
experimental criteria requiring the van't Hoff to calorimetric enthalpy ratio 
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AH vii / AH ca \ ~ 1 (the calorimetric criterion), as well as a two-state-like vari- 
ation of the average radius of gyration upon denaturation. We find that both 
local and nonlocal interactions are critical for thermodynamic cooperativity. 
Chain models with either much weakened local conformational propensities 
or much weakened favorable nonlocal interactions are significantly less coop- 
erative than chain models with both strong local propensities and strong fa- 
vorable nonlocal interactions. These findings are compared with results from 
a recently proposed lattice model with a local-nonlocal coupling mechanism; 
their relationship with experimental measurements of protein cooperativity 
and chain compactness is discussed. 
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1 Introduction 



How a globular protein can fold reliably into a particular three dimensional 
conformation in vitro, without the participation of molecular chaperones, is 
a central puzzle in biophysics. If we wish not only to predict the folded state 
of a protein, but also to understand the folding phenomenon in terms of 
physical processes, we need to use physics-based methods: we run computer 
simulations of self-contained polymer models [Tj that attempt to mimic the 
behavior of real protein molecules. A complete quantum mechanical simula- 
tion, which would include the solvent molecules in addition to all the atoms 
in the protein molecule, is not yet possible. But in attempting to design sim- 
plified models, we face the problem of how to simplify: which characteristics 
are essential and which can be neglected? What effective energy functions 
does this imply for the simplified system? 

Part of this general question is addressed in this article. Folding experi- 
ments on small globular proteins have long shown evidence of thermodynamic 
and kinetic cooperativity 0E|, which indicates a phenomenon similar to a 
first order phase transition between native and denatured states. As our 
group has argued recently [3JIE1IE10IE] 5 this observation can be exploited to 
constrain the set of possible simplified models and interaction schemes: for 
a particular simplified model to be a quantitatively accurate representation 
of protein thermodynamics and kinetics, it is essential that, when appropri- 
ately applied to a small globular protein, it can produce the experimentally 
observed generic cooperative behavior. 

Such constraints turn out to be rather stringent. It is nontrivial to con- 
struct model interaction schemes that can produce proteinlike cooperativi- 
ties lUEJEHZIIHl • A case in point is a class of common Go-like models (HUH]- 
Their potential functions are native-centric, in that they are explicitly biased 
to favor a given native structure. Go-like modeling of proteins has provided 
important physical insights [Tfl | fTT | IT2 t lT3] . These include an increasing num- 
ber of elegant elucidations of functional protein dynamics under native con- 
ditions [IHIinillEllIZllIBl • As for global folding and unfolding of proteins (in 
contrast to their near-native dynamics), a detailed discussion of the merits 
and limitations of Go-like approaches can be found in Ref. |Hj. Notably, 
common Go-like models do not appear capable of producing simple two- 
state folding/unfolding kinetics. Instead, their chevron plots exhibit severe 
rollovers, which are typical of the class of folding kinetics that is customarily 
referred to as non-two-state |H|Hn). Nonetheless, common three-dimensional 
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Go-like protein models seem sufficient to produce apparent two-state ther- 
modynamic behavior [SHE], although their two-dimensional counterparts fail 
to do so |3]. 

In the present investigation, we limit our scope to thermodynamic coop- 
erativity. Specifically, we aim to explore how protein thermodynamic two- 
state-like behavior is affected by the relative strengths of local interactions 
(between residues close together along the chain sequence) and nonlocal inter- 
actions (between residues far apart along the chain sequence). The respective 
roles of local and nonlocal interactions are an issue of long standing interest 
in the study of protein energetics |20 [l211l22t l23 [l24ll25] . and the effect of anal- 
ogous interactions on the phase diagram of lattice polymers has also been 
investigated [2E1E3EH1I2H]- Here the issue is addressed by varying the po- 
tential function in a series of coarse-grained Go-like models, which represent 
the protein as a string of C a positions in continuum space and which are 
simulated using Langevin dynamics. In view of the limitations of common 
Go models (Hj, the present study should be viewed as a first step in tackling 
the issue of local vs. nonlocal interactions in cooperative continuum protein 
models. To assess the robustness of our conclusions, results from the con- 
tinuum Langevin models are also compared with results from lattice model 
simulations. 

We begin in section 2 by providing details of the models. An outline of 
the thermodynamics involved in interpreting the simulations is given in sec- 
tion 3. Our findings are presented in section 4, and we conclude in section 5 
with a discussion of the implications of our results. 



2 Models and simulation details 
2.1 Continuum models 

For the present continuum Go-like models we use a representation, introduced 
by Clementi et al. [10], of the 64-residue truncated form of chymotrypsin 
inhibitor 2 (CI2). The native contact set corresponds to NCS2 in Ref. [E]. 

We use an energy function that is similar to one used previously [8, 10. 30j. 
The potential energy function V, from which the conformational force is 
derived, is given by 

V Stretching ^bending ^torsion ^native ^nonnative; (1) 
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where 

JV-l 



Stretching — E ~ ^o) (2) 



i=l 

contains a summation over the virtual bonds between pairs of residues, 

N-2 

Sending = E £e (^ ~ ^o) 2 (3) 
i=l 

involves a summation over the virtual bond angles between triplets of residues, 
and 

N-3 

Korsion = E I 1 - COS (/ - ^o)] + £ f [ l ~ COS 3 " <&)] } ^ 

i=l 

represents the virtual torsional potential between quadruplets of residues. 
The latter contains a term with a single minimum as well as the traditional 
three-minimum term [31]. [We note that there is an apparent typographical 
error in the corresponding Vision in [10], which effectively lists these terms 
as 1 + cos ((ft 1 — O ) and 1 + cos 3 (0* — But such terms would fold the 
chain into the mirror image of the PDB structure.] Stretching, Sending arid 
Sorsion together account for the local interactions (between residues that are 
separated by no more than three places along the chain), which include local 
conformational propensities for the native structure. The local interactions 
are expressed in this way because it biases the local geometry of the chain. 
The fourth term, 



r.,... 

l»-j|>4 



native / , 




(5) 



sums over the pairwise interactions between residues that are regarded as 
being in contact in the native structure; this accounts for the nonlocal in- 
teractions (between residues that are separated by four or more places along 
the chain). Finally, 

/r re \ 12 

Sionnative / , £ ( JJ~ ) (6) 

contains repulsive pairwise interactions between other pairs of residues, in 
order to ensure the self-avoidance of the chain. 
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The chain contains N residues. I 1 is the length of virtual bond i, 9 l is a 
bond angle, (p 1 is a dihedral angle and r y is the distance between two residues 
i and j. The corresponding values in the native structure are l l , 9 l , (f) l and 
Tq . The range r rep of the repulsive potential between pairs of residues that 
are not bonded and do not interact via a native contact interaction is set to 
4 A. Length is expressed in units of A, and energy in units of e, the energy 
parameter of the nonnative repulsive interaction, so that e itself is unity, ki, 
Be, ejj and £ native are also parameters of the potential energy function. 

hi is fixed at lOOe/A, but the other parameters can be varied; Eq = 20e^ 

and = 0-5e^ are defined in terms of which can be varied to test 
the effect of changing the strength of the local interactions, while £ na tive can 
also be varied (see below) in order to test the effect of changing the strength 
of the non-local interactions. The energy of the system is thus controlled by 
three parameters e, £ na tive and e} . All interaction parameters are taken to 
be temperature independent in the present study. 

Apart from the variable parameters, this energy function differs from the 
similar function used in 0103121] in two further important ways, as follows. 



(1) For r u /ro < y5/6, we set Native = £, while for r 13 /t % q > y5/6, we 
set ^native = £ a - Then we can vary the native interaction parameter e OJ in 
order to test the effect of changing the strength of the non-bonded attractive 
interactions between residues, while the short-range repulsive part of Votive 
maintains the self-avoidance of the chain. (2) The value q, which is the 
smallest number of places along the chain by which two residues can be 
separated if they are to interact by ^normative, can be set either to q = 4 (in 
order to eliminate any double counting of local interactions, in situations 
where is not being varied) or to q = 2 (in order to allow to decrease 
without compromising the self-avoidance of the chain). 
The equation of motion of each residue is 

dvHt) 

where m is the mass of a residue (set to unity), 7 is the coefficient of friction, 
t is time, and v % (t) s F* oni (t) and rf{t) represent each of the three components 
of the velocity, conformational force and random force, respectively [S2J - The 
random force is given by 



•fit) = tH^, (8) 
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where St is the integration time step and £ l is a random variable taken from 
a Gaussian distribution with zero mean and unit variance. The most appro- 
priate time scale can be estimated [H2] by r = ^m Q aQ/eo, where mo, ao and 
e are the mass, length and energy scales, respectively. We set m = m = 1, 
Eq — e — 1 and a = 4 A (the latter is approximately the length of a virtual 
bond between two residues and is also the range r rep of the repulsive interac- 
tion), and so r = 4. We define the integration time step St = 0.005r and the 
coefficient of friction 7 = 0.05r -1 in terms of this time scale. The velocity- 
verlet algorithm [BUSSES] is used to integrate the equations of motion. 

2.2 Lattice models 

The lattice models considered here are 27mers with a maximally compact na- 
tive (ground-state) conformation. Details of the models have been described 
elsewhere |34|I35|. We compare three native-centric interaction scenarios 
which have varying degrees, and different mechanisms, of thermodynamic 
cooperativity. As an example of a particular native conformation to which 
these three scenarios can be applied, we choose the one in Ref. [HE] with rel- 
ative contact order 0.410. In scenario (i), which corresponds to the common 
Go model, the native contact interactions are pairwise additive. In scenario 
(ii), we add an extra favorable energy E gs for the native structure as a whole 
(as defined by equation 5 in Ref. f34j). Scenario (iii) introduces, in place 
of the extra favorable energy, a coupling between the strength of the con- 
tact interaction and the local geometry: two residues which are in contact in 
the native state will interact strongly only when the local geometries of the 
protein chain around the residues are the same as those in the native state, 
as described in Ref. [3E1- We characterize this mechanism as local-nonlocal 
coupling or "a cooperative interplay between favorable nonlocal interactions 
and local conformational preferences" [22]. In this scenario, the strength of 
the native contact interaction is reduced by an attenuation factor a when the 
local geometry is nonnative. The common (uncoupled) Go model is equiva- 
lent to a = 1 (no attenuation), while a = implies complete coupling; a = 
is used here. Under scenarios (i) and (iii), the native state has an energy of 
—28 units, while the extra favorable native energy in scenario (ii) changes 
the energy of the native state to —42 units. Standard Monte Carlo methods 
are used for conformational sampling |34|I35| . The permitted chain moves are 
end flips, corner flips, crankshafts and rigid rotations. Each attempted move 
is counted as one simulation time step, irrespective of whether the move is 
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accepted by the Metropolis criterion. 



3 Thermodynamics 

All simulations are performed at constant temperature, with no explicit con- 
sideration of pressure. This is because the focus of the present study is 
protein behavior under atmospheric pressure, and the contribution of a PV 
term to protein energetics is small under these conditions [3]. Therefore, for 
our present purposes, we can consider the Helmholtz and Gibbs free energies 
to be equivalent. 

3.1 Calculation of the heat capacity 

The specific heat capacity Cy(T) of the model protein is given by the stan- 
dard relation 



where k^F is the Boltzmann constant multiplied by the absolute tempera- 
ture, (X(T)} denotes the Boltzmann average of quantity X at temperature 
T, and the total energy E is the sum V+E K of potential and kinetic energies. 
We compute the averages by standard histogram sampling techniques |8|fT0]. 

In lattice studies, the kinetic energy Ek is not treated. Therefore, E in 
Eq. has traditionally been taken, in protein modeling, to be the potential 
energy term V . This procedure has often been extended to continuum model 
studies, in Ref. [Hj for example, although E K is accessible and well-defined 
in off-lattice models. However, (E 2 ) — (E) 2 ^ (V 2 ) — (V) 2 in general. The 
two quantities would be equal if Ek were a constant, but that would be 
unphysical. In this study, we have calculated CyiT) using Eq. (jOJ both with 
E = V + E K and with the substitution E — > V. The results are not identical: 
an example is given in Fig. 1. 

Fig. 1 shows that the difference between heat capacity values obtained 
using the two methods is small around the transition midpoint T m . This 
is because any energy added to the system during the unfolding transition 
contributes mostly to the potential rather than to the kinetic energy. The dif- 
ference is less negligible for the "shoulders" on either side of the heat capacity 
peak. At very low temperatures, including the kinetic energy contribution 





8 



can lead to a smaller heat capacity, because the molecule at this tempera- 
ture is in a relatively fixed state: nearly all of the kinetic energy is accounted 
for by the oscillation of pairs of residues about the minima of their mutual 
(bonded or non-bonded) interaction energies. The potential energy and the 
kinetic energy associated with these oscillations both fluctuate, but their sum 
fluctuates much less, and so the fluctuations in total energy are smaller than 
the fluctuations in potential energy, with the result that the calculated heat 
capacity is smaller when Ek is taken into account. Overall, Fig. 1 indicates 
that while the difference between the heat capacities calculated using the 
two different methods is not negligible, it is not drastic. Probably this is 
because Ek, while not invariant, fluctuates much less than V. For this rea- 
son, we do not expect conclusions drawn from previous calculations of heat 
capacities [E], which used V, to be changed greatly by calculations using 
E. Nonetheless, we do expect a proper account of the kinetic contributions 
to protein heat capacities to be important in addressing the contribution of 
bond vector motions to the heat capacity [36] . 

Since the PV term is neglected in the present formulation, Cy(T) is 
effectively equal to Cp(T), which is generally measured by calorimetry (and 
which can be expressed in a form similar to Eq. (jHJ), but with the enthalpy 
H taking the place of the energy E |4,5j). Therefore, we may refer to the 
quantity computed using Eq. (JUJ simply as heat capacity. All subsequent heat 
capacity curves shown in this article for the continuum models are obtained 
using the total energy E = V + Ek- 



3.2 The free energy 

The Helmholtz free energy of the model system is F(T) = — k B Tln Z(T), 
where Z(T) is the partition function at temperature T. It follows that, in 
the vicinity of the simulation temperature T sim , the Helmholtz free energy of 
the model system at temperature T, relative to its value at the simulation 
temperature, may be approximated using the formula: 



k B T k B T k B T sim IV 1 



1 



- k B T s \ m k B T ^ 
(10) 

where the sum is performed over sets of microstates in different energy ranges 
Ei. p(Ei] T sim ) is the probability density at the simulation temperature, and 
is estimated directly from the Langevin dynamics simulations. 
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The inset of Fig. 1 provides an example of AF(T), showing that the gra- 
dient of the free energy with respect to T changes rather abruptly around 
the transition temperature T m (vertical dotted lines). The transition tem- 
perature T m corresponds to the temperature at the peak of the heat capacity 
curve, which was denoted by T max in Ref. Apparently, below T m , the 
system spends most time in states in the vicinity of the bottom of the native 
basin, and so the changes in F with respect to temperature are dominated 
by the behavior of these states. However, as the temperature increases past 
T m , the system, and therefore the rate of change of F(T), starts to be dom- 
inated by states near the bottom of the denatured basin. As a result, the 
gradient of F(T) changes rather suddenly at T m . While the F(T) gradient 
could never be discontinuous because the model system is finite, the kink at 
T m does indicate that the transition is two-state-like, and therefore that it is 
similar to a first order phase transition. 

3.3 Thermodynamic cooperativity 

The presence of a peak in the heat capacity at a transition temperature 
T m , as in Fig. 1, indicates that the folding/unfolding transition possesses 
a degree of thermodynamic cooperativity. As our group has argued, the 
degree of thermodynamic cooperativity in protein models can be quantified 
by the ratio H2 = AH v n/ AH ca \ of the van't Hoff enthalpy Aif v H to the 
calorimetric enthalpy AH ca \ of the transition. This ratio is closely related to 
that determined experimentally by differential scanning calorimetry |37| . In 
model studies, the calorimetric enthalpy AH ca i may be determined from an 
integral of the heat capacity across the transition region, 



while the van't Hoff enthalpy is equal to twice the maximum standard devi- 
ation of the enthalpy distribution at the transition midpoint, 



where Cp ax is the peak value of the heat capacity. We have followed standard 
usage in this section by expressing k 2 in terms of H and Cp. However, as 
mentioned in the previous section, simulations produce values for E and Cy, 
which for the present application are essentially equivalent to H and Cp. 




(11) 




(12) 
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As has been pointed out [E], comparison of simulation heat capacity scans 
to experiment is often complicated by the fact that the heat capacity tails 
which we observe in simulations would, if they occurred in a real system, be 
swamped by the solvent contribution and ignored by the common procedure 
of using empirical baseline subtraction to calculate AH v yi/ AH ca \. In other 
words, tail contributions that arise from conformational transitions may be 
masked by solvent contributions in real data analysis [SllZj- Therefore, for 
completeness, we also perform empirical baseline subtractions on our simu- 
lated heat capacity scans, producing a revised ratio /4 (defined in Ref. [H]) 
to facilitate comparison with experiment. 

3.4 Radius of gyration 

The radius of gyration R g of a particular conformation of the protein is an 
indicator of its compactness. It is defined by 



where N is the number of residues, is the position of the ith residue, 
and (r) is the average position (centroid) of the given conformation. The 



obtained by standard histogram techniques. Two-state-like behavior requires 
a steplike sigmoidal change in (R g ) upon denaturation at T m , with little 
postdenaturational expansion of the chain [5]. 

4 Results and discussion 

To study the effect of local vs. nonlocal interactions, we first vary the strength 
of the local interactions while keeping the strength of the nonlocal in- 
teractions fixed in the continuum CI2 construct (Figs. 2, 3). The heat ca- 
pacity scans, for four scenarios (four models) with different values for el\ 
are shown in Fig. 2. They all exhibit a fairly sharp peak except for the 
model with = 0.25. The heat capacity peak signifies substantial heat 
absorption within a narrow temperature range at the folding/unfolding tran- 
sition. The absorbed energy propels the chain from its low-energy folded 
conformations (native ensemble) to its high-energy unfolded conformations 



1 N 
iV i=i 



(13) 




over a given conformational ensemble is 
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(denatured ensemble). However, the mere existence of a relatively sharp peak 
in the heat capacity function does not necessarily mean that the transition 
is as cooperative as those observed in small single-domain proteins. Coil- 
globule transitions in homopolymers are not two-state-like, but their calori- 
metric heat capacity scans can have very sharp peaks [38j. A more quantita- 
tive measure of thermodynamic cooperativity is the traditional calorimetric 
two-state criterion (see above), which has emerged recently as a powerful 
modeling tool HIElEll3IHlEnilinillllll2j ■ Ratios of van't Hoff to calorimetric 
enthalpy were calculated for the four models as described above; the ranges 
of the temperature integrations used in the determination of AH ca \ [Eq. (fTTJ) ] 
were taken to be equal to the ranges shown in Fig. 2. 

The inset of Fig. 2 shows that the AH v yi/ AH ca \ ratio (diamonds) of these 
models is only weakly dependent on over an extended range of val- 
ues, but that the ratio is significantly smaller when the local interactions 
are substantially weaker, at e$ = 0.25, than the nonlocal interactions. As 
discussed above, quantitative comparisons between simulated and experimen- 
tal Aif v n/AiJ ca i values require the introduction of model calorimetric base- 
lines similar to those employed in the interpretation of experimental data. 
Traditionally, experimental baselines are designed to remove solvation con- 
tributions (temperature-dependent effective interactions), in order to extract 
the heat capacity effects associated with the folding/unfolding transition it- 
self [37J. The present models do not contain temperature-dependent inter- 
actions. Therefore, the heat capacity contributions eliminated by the model 
calorimetric baselines in Fig. 2 can only originate from vibrational motions 
and conformational transitions. Increasingly, it is being recognized [5|l36|l^3] 
that similar heat capacity contributions from bond vector motions and more 
collective conformational transitions might also be "hidden" below traditional 
baselines constructed for analyzing experimental calorimetric data, although 
the magnitude of such contributions needs to be elucidated. The three mod- 
els in Fig. 2 that are relatively more cooperative (with higher k 2 values) all 
have modified Aif vII / AH cal values (k^ \ circles in the inset), after empirical 
baseline subtractions, that are very close to unity. (We note that the recent 
determination of AiJ vH /Aif ca i values in an all-atom Go model [l2j involved 
baseline subtractions as well: c.f. Fig. 8 of Ref . |42j.) 

However, a protein chain model's ability to attain a near-unity AH v u/ AH ca \ 
ratio after baseline subtractions does not by itself imply that its thermody- 
namic behavior is similar to that of real, small single-domain proteins (HIE]- 
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This is because the heat capacity contributions discarded by certain baselines 
can actually be symptoms of significant deviations from two-state-like behav- 
ior. It has been recognized [2] that, to clarify this situation, we can use the 
behavior of the average radius of gyration (R g ) of a protein chain model as an 
additional evaluation criterion for the model's thermodynamic cooperativity. 
Small angle X-ray scattering (SAXS) experiments have demonstrated that 
the average radius of gyration (R g ) of several small single-domain proteins 
behaves in an apparently two-state manner jUJHHHIH], showing very little 
postdenaturational (T > T m ) expansion of the chain outside the transition 
regime that corresponds to the region of the heat capacity peak. We require 
chain models of small single-domain proteins to exhibit similar behavior [H]. 
Now, to further assess the four models in Fig. 2 with different local interac- 
tion strengths, we calculate their average radii of gyration as a function of 
temperature (Fig. 3). To ensure adequate sampling, (R g ) for each model is 
obtained from three different simulation temperatures; the results are thus 
displayed as three discontinuous curves. Despite some minor discrepancies 
(owing to sampling uncertainties) between parts of the (Rg) function deduced 
from different simulation temperatures for the = 0.25 case, the general 
trend in Fig. 3 is very clear. Models with weaker local interactions are less 
cooperative in that their (R g ) curves show more postdenaturational increase 
than do those of models having stronger local interactions. For instance, 
the (Rg) of the ef = 0.25 model increases by « 3.0 A between T « 0.82 
(the end of the transition region) and T ~ 1.11. In contrast, a similar tem- 
perature increase for the = 1.00 model from T m 1.11 (the end of the 
transition region) to T rs 1.40 leads to an increase of only ~ 1.6 A in (R g ). 
These observations indicate that two-state-like thermodynamic cooperativ- 
ity cannot be achieved if the local conformational propensities of a protein 
are much weaker than the favorable nonlocal interactions. This confirms a 
similar conclusion which was derived recently from a more limited study of 
a "contact dominant model" [E|. 

We next extend our analysis by applying the same computational proce- 
dure to varying the strength e a of the favorable nonlocal interactions while 
keeping the strength of the local interactions fixed. Consistent with the sem- 
inal study of Go and Taketomi (2D], Figs. 4 and 5 show that variations in 
nonlocal e a have a more prominent effect on thermodynamic cooperativity 



models in Fig. 2 with e\ > 0.5 are similar, the peak heat capacity values for 




While the peak heat capacity values for the three 
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the three models in Fig. 4 with e a > 0.5 show a significant monotonic increase 
with e a . In addition, for the e a = 0.5 model in Fig. 4, the difference between 
unity and the AH vi i/ AH cal ratio after baseline subtraction is not negligible 
(«2 = 0.91). Despite these differences, the trends in Figs. 4, 5 are in large 
measure similar to those in Figs. 2, 3. In particular, Fig. 4 shows that the 
model with e a = 0.25, like the ef = 0.25 case in Fig. 2, has a significantly 
lower AH v u/ AH ca \ ratio than the other three models considered in the same 
figure. The (R g ) data in Fig. 5 shows that thermodynamic cooperativity 
increases with e a , as manifested in a smaller amount of postdenaturational 
conformational expansion with increasing e a ; this is comparable to the effect 
of increasing in Fig. 3. Taken together, the results in Figs. 2-5 suggest 
that a high degree of thermodynamic cooperativity, similar to that in real, 
small single-domain proteins, requires both strong local and strong nonlocal 
interactions. Apparently, a high degree of thermodynamic cooperativity is 
incompatible with either a much weakened local conformational preference 
relative to the favorable nonlocal interactions {er}' <C e a ) or much weakened 
favorable nonlocal interactions relative to the local conformational preference 

(e.<4 1) )« 

Although three-dimensional Go-like models with strong local and nonlo- 
cal interactions appear to satisfy the thermodynamic criterion of calorimetric 
two-state cooperativity, it has recently been noted that they are unable to 
produce simple two-state folding/unfolding kinetics [7|l8] . This is because the 
thermodynamic cooperativity of these models is not sufficiently high. As a re- 
sult, and in spite of the native-centric nature of the common pairwise additive 
Go-like interactions, kinetic trapping becomes significant under strongly na- 
tive conditions, leading to folding rate slow-downs and chevron rollovers [19j. 
More recent lattice model investigations indicate that simple two-state fold- 
ing/unfolding kinetics require a high degree of thermodynamic cooperativity 
that may be characterized as "near-Levinthal" [HI], necessitating many-body 
interactions beyond those postulated by the common Go model [3i|l35]. 

In view of this recent development, and to facilitate the construction and 
investigation of continuum models that incorporate these new ideas, it is in- 
structive to compare in more detail the thermodynamic behavior of the com- 
mon lattice Go construct (with only pairwise additive contact energies) with 
the behavior of models which have many-body interactions and enhanced 
cooperativity. We also wish to investigate whether results obtained from 
lattice models supply additional support for the conclusions which we have 
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derived from our continuum model results. To this end, Figs. 6-8 compare 
three 27mer lattice models. 

Because of their intrinsic restrictions on conformational possibilities, it 
is more straightforward to construct cooperative lattice models than to con- 
struct off-lattice continuum models that are similarly cooperative. Recently, 
using evidence from kinetic simulations of chevron plots, our group has pro- 
posed that a AH vi {/ AH cal ratio of k 2 > 0.9 before baseline subtractions [as 
for models (ii) and (iii) in Fig. 6a] is likely to be required in order for a lattice 
protein chain model to produce chevron plots with linear regimes similar in 
extent to those observed for real, small single-domain proteins P3j. However, 
this numerical criterion is not readily generalizable to off-lattice continuum 
models. This is because the heat capacity effects of bond vibrations and 
kinetic energy have to be taken into account in continuum models, whereas 
these effects are absent in lattice models. Thus, in the characterization of 
a model's thermodynamic cooperativity, more detailed information concern- 
ing, for example, the behavior of the average radius of gyration, has to be 
relied upon more heavily for continuous models (see above) than for lattice 
models. 

For the three lattice models studied here, the (R g ) plots in Fig. 6b show 
that, while the postdenaturational conformational expansion of the common 
lattice Go construct (solid curve in Fig. 6b) is considerably milder than that 
of its continuum counterpart (solid curves in Figs. 3, 5), the more cooper- 
ative lattice models with many-body interactions exhibit much less (dotted 
curve in Fig. 6b) or nearly non-existent (dashed curve in Fig. 6b) postdenat- 
urational conformational expansion. The fluctuations in E and R g near the 
transition midpoint, shown in Figs. 7 and 8, indicate further that the tran- 
sitions between the native and denatured ensembles are sharper and more 
two-state-like for the more cooperative models (ii) and (iii) with many-body 
interactions [parts (b) and (c) of Figs. 7, 8] than for the common Go model 
[parts (a) of Figs. 7, 8]. The corresponding fluctuations in the fractional 
number of native contacts Q [34, 35| (data not shown) were also found to 
exhibit a trend very similar to that of the energy fluctuations in Fig. 7. 

The lattice model results, shown in Figs. 6-8, are compatible with the con- 
clusion, reached above on the basis of continuum model results, that both 
local and nonlocal interactions are important for thermodynamic coopera- 
tivity. The common lattice Go model of scenario (i) includes only nonlocal 
interactions, analogous to the interactions encoded by Kativc in the off-lattice 
model. Scenario (iii), which takes account also of the local geometry of the 
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chain, displays greater cooperativity than scenario (i). 

Higher resolution data such as that in Figs. 7 and 8 opens up future av- 
enues for the assessment of different mechanisms of cooperativity using com- 
parisons between model predictions and experimental measurements of, for 
example, conformational sizes and fluctuations - It is noteworthy 

that in the model (ii) scenario, with an extra favorable energy for the native 
structure as a whole, the native ensemble does not exhibit much energetic or 
conformational variation (horizontal line segments at low E and low R g val- 
ues in Figs. 7b and 8b). On the other hand, in the model (iii) scenario with 
local-nonlocal coupling, there is considerable variation in the native ensemble 
(c.f. low E and low R g fluctuations in Figs. 7c and 8c). Yet the variation in 
the denatured ensemble is smaller in model (iii) than in model (ii) (c.f. high 
E and high R g fluctuations in parts (b) and (c) of Figs. 7, 8), resulting in 
a more two-state-like average (R g ) transition for model (iii) than for model 
(ii), manifested in a near-immediate postdenaturational (T > T m ) saturation 
of the dashed curve for model (iii) in Fig. 6b compared to a more gradual 
postdenaturational saturation of the dotted curve for model (ii) in the same 
figure. All these differences in conformational properties are in principle de- 
tectable through experiments on real proteins. Hence, future experimental 
efforts along the lines suggested here would help to verify or falsify different 
proposed scenarios and interaction mechanisms [34ll35U50|l5T] in the endeavor 
to decipher the physical origins of cooperativity in real proteins. 

5 Conclusions 

The present study suggests strongly that both local conformational prefer- 
ences and favorable nonlocal interactions are crucial to protein thermody- 
namic cooperativity. This result points to a useful constraint on simplified 
models of protein molecules: they should take account both of local and of 
nonlocal interactions. 

As emphasized above, the scope of the present study is limited. Only 
native-centric interaction schemes are considered; and here we have elected 
only one particular physically plausible way to classify energy contributions 
into "local" and "nonlocal" terms in the continuum models. In addition, 
by using a native-centric interaction scheme both for local and for nonlocal 
interactions, we have avoided the possibility of energetic frustration, which 
might be significant if a more realistic interaction scheme were used [H21 
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153] . To further elucidate the answers to the questions we have posed, much 
remains to be investigated. 

Nonetheless, our results show clearly that a high degree of thermodynamic 
cooperativity is compatible neither with a much weakened local interaction 
nor with a much weakened nonlocal interaction, indicating that both local 
and nonlocal interactions are important components in protein energetics 
[Hi] , This finding is consistent with the notion that a cooperative interplay 
between local and nonlocal interactions |HE30Elll2E| is a critical ingredient 
underlying the apparent simple two-state cooperativity of real, small single- 
domain proteins. 

With regard to Go-like native-centric modeling (see, e.g., discussion in 
Refs. [H1III3HHEI3EE1EI!), we 0Dserve that significant differences in model 
predictions can result from different Go-like interaction schemes, even though 
all of the schemes are designed to bias the chain towards the same native 
structure. This underscores our point that requiring a consistent account 
of cooperativity can be a more productive approach to protein modeling 
than simply designing a model heteropolymer to fold to a target structure 
[B|. In this context, the present coarse-grained representations constitute 
only a first step in the understanding of protein cooperativity. Ultimately, 
atomistic origins of local and nonlocal interactions such as sidechain packing 
[58ll59U60|l6T] must be taken into account in an effort to provide the necessary 
physical underpinning for the cooperative mechanisms proposed here. 
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Figure Captions 

FIGURE 1 

Heat capacity as a function of temperature, for one particular set of in- 
teraction parameters. Solid curve: heat capacity calculated using the total 
energy; dashed curve: heat capacity calculated using only the potential en- 
ergy. Parameters are = 1.00, e a = 1.00, q = 4; simulation temperature 
^sim — 1-02. Inset: free energy as a function of temperature, for the same 
model, showing a sharp change in gradient around T m . The vertical dotted 
lines in the figure and the inset mark the transition midpoint temperature T m . 

FIGURE 2 

Heat capacity as a function of temperature, for varying local interac- 
tion energy el\ Other parameters e a = 1.00 and q = 2 are fixed. (From 

left to right) dotted curve: = 0.25, T m = 0.74; short dashed curve: 

e W = 0.50, T m = 0.84; long dashed curve: = 0.75, T m = 0.94; solid 

curve: e^P = 1.00, T m = 1.03. These scans are obtained by histogram 
techniques from simulations performed at T sim = 0.73, 0.84, 0.94, and 1.03 
respectively. The AH v n/ AH ca \ cooperativity coefficients k 2 without baseline 
subtractions are 0.33, 0.43, 0.44, and 0.44 respectively. Modified cooperativ- 
ity coefficients after subtraction of the baselines (indicated by thin lines 
in the figure) for e£ } = 0.50, 0.75, and 1.00 are 0.97, 0.98, and 0.99 respec- 
tively. No value for was calculated for = 0.25 because the shape of 
its heat capacity curve does not suggest any clear choice of baselines that are 
intuitively more reasonable than others. The inset shows n 2 (diamonds) and 
(circles) as functions of 

FIGURE 3 

Average radius of gyration as a function of temperature, for varying local 
interaction energy el , other parameters e a = 1-00 and q = 2 are fixed, as 

in Fig. 2. The correspondence between line styles and values is identical 

to that in Fig. 2. For each value of \ simulations were performed at three 
different values of T sim to ensure adequate sampling across the entire tem- 
perature range shown. (From left to right) for = 0.25 (dotted curves), 

T sim = 0.73, 0.93, 1.13; for e£ } = 0.50 (short dashed curves), T sim = 0.84, 

1.04, 1.24; for e { P = 0.75 (long dashed curves), T sim = 0.94, 1.14, 1.34; and 
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for 4 1} = 1-00 (solid curves), T sim = 1.03, 1.23, and 1.43. 
FIGURE 4 

Heat capacity as a function of temperature, for varying nonlocal interac- 
tion energy e a . Other parameters = 1.00 and q = 4 are fixed. (From left 
to right) dotted curve; £ a — 0.25, T m — 0.43; short dashed curve; £ a — 0.50, 
T m = 0.65; long dashed curve: e a = 0.75, T m = 0.84; solid curve: e a = 1.00, 
T m = 1.02. These scans are obtained by histogram techniques from sim- 
ulations performed at T sim = 0.42, 0.64, 0.84, and 1.02 respectively. The 
AH v u/ AH ca \ cooperativity coefficients k 2 without baseline subtractions are 
0.28, 0.40, 0.44, and 0.46 respectively. Modified cooperativity coefficients 
after subtraction of the baselines (indicated by thin lines in the figure) 
for e a = 0.50, 0.75, and 1.00 are 0.91, 1.00, and 0.99 respectively. No value 
for was calculated for e a = 0.25 for the same reason that no was 
provided for = 0.25 in Fig. 2. The inset shows k 2 (diamonds) and nf 1 
(circles) as functions of e a . 

FIGURE 5 

Average radius of gyration as a function of temperature, for varying non- 
local interaction energy e a . Other parameters = 1.00 and q = 4 are 
fixed, as in Fig. 4. The correspondence between line styles and e a values is 
identical to that in Fig. 4. For each value of e a , simulations were performed 
at three different values of T sim to ensure adequate sampling across the whole 
range shown. (From left to right) for e a = 0.25 (dotted curves), T sim = 0.42, 
0.62, 0.82; for e a = 0.50 (short dashed curves), T sim = 0.64, 0.84, 1.04; for 
e a = 0.75 (long dashed curves), T sim = 0.84, 1.04, 1.24; and for e a = 1.00 
(solid curves), T sim = 1.02, 1.22, and 1.42. 

FIGURE 6 

Thermodynamic cooperativities of the three representative 27mer lattice 
models described in the text. Model definitions and simulation details are 
given in Refs. [29,30]. Heat capacity (a) and average radius of gyration (b) 
are determined by histogram techniques based upon Monte Carlo sampling 
performed at T sim = T m . (a) Heat capacity of (i) the common Go model 
with pairwise additive contact energy (solid curve, left); (ii) the model that 
assigns an extra favorable energy to the native structure as a whole (dotted 
curve, right); and (iii) the model with local-nonlocal coupling (dashed curve, 
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middle). The transition temperatures for the models are T m = 0.701 (i), 
1.13 (ii), and 0.755 (iii). Their AH v u/ AH ca \ without baseline subtractions 
are k 2 = 0.86, 0.98, and 0.99, and the corresponding ratios after subtracting 
the baselines shown are = 1.00, 1.00, and 1.00 respectively, (b) Average 
radius of gyration as a function of model temperature for the three models 
[represented by the same line styles as in (a)]. 

FIGURE 7 

Representative trajectories of the three models in Fig. 6 at their respec- 
tive transition temperatures. Variations in the potential energy of models 
(i)-(iii) are shown in (a)-(c) respectively. 

FIGURE 8 

Same as Fig. 7, except that variations in the radius of gyration are shown 
here. 
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